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In the context of equation-free computation, we devise and implement a procedure for using 
short-time direct simulations of a KPZ type equation to calculate the self-similar solution for its en- 
semble averaged correlation function. The method involves "lifting" from candidate pair-correlation 
functions to consistent realization ensembles, short bursts of KPZ-type evolution, and appropriate 
rescaling of the resulting averaged pair correlation functions. Both the self-similar shapes and their 
similarity exponents are obtained at a computational cost significantly reduced to that required to 
reach saturation in such systems. 

I. INTRODUCTION 

Often we are faced with the study of systems described by a given fine scale, microscopic dynamics, for which 
we would like to obtain coarse-grained, macroscopic information. Such information can include stationary-states, 
instabilities, and bifurcations. When continuum equations describing the coarse-grained dynamics are available in 
closed form, traditional numerical analysis tools can be used to obtain this type of information efficiently. Recently, 
much work has been devoted to developing tools for addressing such questions in the absence of an explicit coarse- 
grained description (see 0,0.13,111 an< ^ references therein). These equation-free methods offer the hope of significant 
savings in storage and run-time costs over direct numerical simulations of the underlying microscopic dynamics. 
They also allow the study of questions inaccessible to direct simulation, for example the characterization of unstable 
stationary states. These techniques have been applied to the coarse-grained study of systems described by microscopic 
evolution rules (kinetic Monte-Carlo, Brownian and molecular dynamics, Lattice-Boltzmann etc., see the references 
in 3). 

These tools can also be applied to problems whose solutions are macroscopically characterized by a symmetry group, 
such as translational invariance (giving rise to traveling wave solutions). More recently, the techniques have been 
extended to problems whose macroscopic dynamics exhibits scaling; among these are molecular diffusion |{|, models 
of self-similar transport of random Brownian particles [||, core collapse in stellar systems Q etc. Such problems present 
challenges, both conceptual and technical, not necessarily found in the simpler traveling-wave problems. For example, 
since scale-invariance -unlike translational invariance- is never an exact symmetry of the microscopic dynamics, the 
scaling solutions only exist in an asymptotic limit of large system sizes and long time-scales. Furthermore, one must 
correctly identify the number of free parameters characterizing the scaling solution. 

In this context we consider the scaling behavior of an equation in the KPZ class This model is a paradigm for 

a wide class of systems whose correlations exhibit asymptotic scaling. Furthermore, as opposed to the systems studied 
through equation-free methods to date, the scaling behavior is not present in the first moment of the evolving field, 
but rather in its correlations; in particular, we will study the scaling solutions for the two-point correlation function. 
We demonstrate, using the ideas outlined above, how to determine the correlation function self-similar shape and the 
scaling exponents that characterize it. 

The paper is organized as follows: in Section 2 we present the exact form of the equation to be solved, define 
its correlation function and briefly review its known scaling properties. Following that, in Section 3 we discuss an 
iterative method for the determination of the self-similar solution and its exponents. Section 4 presents a matrix-free 
fixed-point approach to the same problem, and we conclude with a discussion in Section 5. 
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II. PRELIMINARIES 



The equation we will analyze herein is a discretized form of a modified KPZ equation (in 1+1 dimensions) for the 
height h of an interface: 

h(x,t) = h"(x,t) + - Xk 't) t)2 ^ +V(x,t), (1) 



l + /j,h'(x,tf 



where the noise rj is (^-correlated in space and time: 



(ry) = 

{rj(x,t)t](x , ,l/)) = SS(x - x')8(t - t'). (2) 

The fi term is present to ensure that the equation does not exhibit the finite-time singularity known to arise for 
sufficiently large A in the unmodified fi = equation 10]. Dimensional analysis shows that // does not change the 
scaling behavior of the system. 

In discretized form, our equation reads 



,t+At 



= h t i +At 



h l+1 2h t + + 4 + ^ h t +i _ h t_ x y + V, 



(3) 



with (rjl) = 0, {rjlrfj ) — S5 tt t>Si,j / At. We work in a periodic box of length L. In the following, we take fx = 1, 
S = 1/12, At = 0.05. It is customary to start simulations (t — 0) from a flat interface h® = 0. It is convenient 
to define the average height h l = i ^2f=Q h\ and the reduced height h\ = h\ — h l . Since a flat interface is stable, 
(hi) = 0, and the basic object of interest is the two-point correlation function: 



G t (d) = -J2(HH +d ) (4) 



or, equivalently, its (discrete) Fourier transform: 

G^ = ^e 27r ^/ L G t (d) (5) 

where K = kL/ (2ir) = 1, . . . , L/2 and k is the momentum. 

Let us briefly review the basic known scaling properties of G l K 0. For intermediately large wavenumber K, 
K*(t) <C K <C L/2, G l K falls like a power-law for increasing K, G l K ~ C/K 2a+1 , where C is independent of t. At 
small K <C i£f*(£), the power-law growth is cut off, and G approaches a finite limiting value as K — > with zero 
derivative. The crossover length scale, l/K*(t), grows with time as t 1 ^, until the saturation time when l/K*(t) ~ L. 
This saturation time thus grows with the system size as L z . At small K K*(t), G saturates at a value which grows 
in time like £( 2a+1 )/ z , again until the saturation time. The values of the exponents a and z in this one-dimensional 
system are known analytically , a = 1/2, z = 3/2. Putting this all together, G has the scaling form 

G l K =t {2a+1)/z f{Kt 1/z ). (6) 
where the function f{x) approaches a constant for small arguments and decays as l/x 2a+1 for large x. 



III. SOLUTION BY DIRECT ITERATION. 



Our goal is to find a self-consistent G l K with the correct scaling properties, i.e. find the function /(.) in Eq. 6 
above. Roughly speaking, if we knew the exact G l K , we could use this knowledge to generate an ensemble of initial h 
fields at some time to, conditioned on this G. We could then evolve each of these initial h's forward in time to some 
tf, and calculate G at this later time. The new G should then simply be a rescaled version of the original G. This 
scaling condition is what we use to determine G. 

There are a number of caveats that need to be expounded at this point. Some of these are technical in nature, 
involving the details of the algorithm used to actually find G. Two, however, involve matters of principle. The first 
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caveat is that the full statistical solution is not uniquely determined solely by the two-point function G. In principle, 
a complete solution implies knowledge of all higher-order correlations as well. However, we posit that the dynamics 
is such that the higher-order correlators relax much more quickly than the two-point function itself. If this is true, 
there will be a short transient during which time the system will reconstruct all the higher-order correlations (will 
"heal") while G itself does not change much. In effect, this is an assumption of separation of time scales between 
different order correlators. Thus, to determine the solution, we compare not G at times to and tf, but at times ti 
and t f , where i/ is an intermediate time chosen so that the higher-order correlators have had time to become slaved 
to G (recover their correct values conditioned on the given G). 

The second caveat is that G exhibits the desired scaling properties only asymptotically. At the smallest scales, the 
underlying lattice ruins scale-invariance. More serious, however, is the fact that for short times and small scales, the 
scaling properties of G are those of the Edwards- Wilkinson, i.e., the A = model Only for sufficiently large 

times is the interface sufficiently rough for the nonlinearity to determine the scaling. Starting simulations from a flat 
interface and evolving for time to brings the interface to some scale of roughness parametrized by to- Our procedure 
will converge on the particular member of the scaling family that has this characteristic roughness scale; needless to 
say, all members of the scaling family (at large enough scales !) can be directly reconstructed. As we will see in more 
detail below, unless the initial time to (alternatively, the initial roughness scale) is large enough, our results will be 
contaminated by the short-time, small scale Edwards- Wilkinson behavior. Only if our working scale is large enough 
for the asymptotically self-similar behavior to have set in, does our process make sense. In practice this means that 
we must test the working roughness scale (alternatively, the working f ) to confirm that both the self-similar shape 
and the exponents have converged. As we increase to, we will of course have to increase L accordingly, so that we 
continue to capture the full self-similar structure of G. 

We now move on to the technical details of the calculation. The primary issue to address is how to solve the 
fixed-point equation embodying the scale invariance condition. We will do this in two different ways. The first is 
by successive substitution. Here we just perform repeated cycles of forward integration, followed by rescaling to the 
original time (roughness). The second is via a fixed point Newton- type procedure, which we describe in the next 
section. 

Direct iteration proceeds as follows. We start by integrating the system forward from a flat interface at t = 
to some to. We then calculate G K . It is straightforward to generate configurations hi conditioned on this G K (to 
"lift" from G to h). All we have to do is to remember that Gk is the expected value of |/i_r-| 2 . Thus, we generate 
a Gaussian random number with zero mean and variance G K , and multiply by a random phase to obtain an hx- 
An inverse Fourier transform gives us our desired hi. We generate some number, (typically 32,000) of such initial 
configurations and integrate each forward in time to tf, which we take to be tf = 2to, and measure G f K . We also 
measure GL at an intermediate time tj = 3to/2. From these measurements we construct the functions G°(K), G 1 (k) 
and G^(k) (see below for details). We now need to rescale G-f(k) back to the original roughness scale (or time to). 
We do this is two steps. We first determine, for each measuring time, a typical small k scale, ki/2 by the condition 

G(fci/2) = G(0)/2. We next rescale wavenumbers at times tj and tf by the factors — ky 2 /k*'L, respectively; 
this "aligns" their rescaled large scale behavior. We next rescale G 1 and G^ by factors f@ = G /, ^(fcs//^'^)/G°(fcs) 
(where ks — \J L/2 is the geometric mean of the smallest and largest wavenumbers) so that their large wavenumber 
behavior coincides with that of G°. For the self-similar shape, the function G? (k/ f£)/fjL should reproduce our 
original function G (K). In practice, it gives us a new starting point for another round of iteration; in our problem 
the self-similar solution is stable, and this procedure rapidly converges, so that in fact the differences are only due to 
fluctuations. If the problem was truly (as opposed to asymptotically) self-similar, upon convergence to the self-similar 

if if 

solutions, the scaling factors and / G could be used to estimate the effective scaling exponents a and z: 

z = ln(t f /t I )/MKl /2 /K{ /2 )=]n(t f /t I )/Hf f K /f I K) 

n i ( HfLlfh 



Mti/flcJ 



- 1 ( ? ) 



Because the self-similarity is only asymptotic, three successive times are in fact needed to estimate the exponents. This 
procedure constitutes the equation-free implementation of dynamic renormalization (see e.g. the classical references 
[HI Edll as well as our template-based approach discussed in 0,E1)- 
The last point we need to cover is how we convert our measured Gk into a function G(fc). We use a relatively 
low-dimensional description of the function G(k); in particular, we fit a cubic spline to log(GK) throughout the whole 
range, with knots approximately equally spaced in log(K); this provides an adequate description with O(10) degrees 
of freedom. 

An example of the results of this procedure are shown in Fig. ^ where one cycle of the iteration is shown. We see 
that indeed G(k) is recovered to quite good accuracy by our integration and rescaling, except for the very largest fc's, 
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FIG. 1: The iteration cycle for A = 5, L = 400, to = 40. a) G(K,t) for t = 40 and t = 80 
together with G(K, 80) after rescaling of if to align its K 1/2 with G(isT, 40). b)G{K,AQ) and G(JC 
after rescaling of K and after rescaling of both k and h to collapse it onto G(K, 40). 



with their respective Ki/z's noted, 



together with G(K, 80) 




10 



15 



Iteration 



FIG. 2: Calculated values of the exponents z, a vs. iteration for A = 5, L = 400, to = 40. 



where the zero-slope condition imposed by the discreteness of the lattice is evident in the original G, but not in its 
rescaled version. 

We present in Fig. [21 a graph of the estimated exponents as a function of iteration number, for some particular 
value of to and L. We see that the iteration is quite stable, albeit with sizable statistical fluctuations. In Fig. EUa-b), 
we present the cumulative average over iterations of the estimated exponents, as a function of iteration number, for 
various sets of to and L. What is remarkable is that the calculation essentially converges immediately, to the precision 
we can measure. Our procedure allows us work at an intermediate scale where we can see simultaneously the rapid 
saturation of the correlation function at short length scales and the growth at the coarse scales, tracking the shift of 
the crossover between these two regimes while this evolution is still relatively fast; approaching the ultimate shape 
while the interface evolves to coarser scales would take a significantly longer computational time (the coarser the scale, 
the longer the time). It is this "shape evolution at constant scale" that underpins the computational savings of the 
method. We see that increasing to (alternatively, the working roughness scale) brings the measured value of z down 
much closer to its asymptotic value of 3/2. Since the value of z for the Edwards- Wilkinson model is 2, we interpret 
this as indicating the contamination of our measurement by the crossover from Edwards- Wilkinson behavior at short 
scales. In Fig. we present for comparison the data for the Edwards- Wilkinson model. Here we immediately obtain 
values very close to the expected z — 2, a — 1/2, as the only violation to scaling comes from the very short scale 
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FIG. 3: Cumulative average of the estimated exponents (a) z; and (b)a as a function of iteration number for A 
pairs of to, L. 
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FIG. 4: Cumulative average of the estimated exponent z as a function of iteration number for the Edwards- Wilkinson model 
(A = 0) and two sets of to, L. 



lattice effects. 



IV. NEWTON-GMRES FIXED POINT SOLUTION 



A direct iteration procedure is, of course, not guaranteed to converge; the self-similar solution may not be stable for 
the system and parameter values of interest. Even if it is stable, one needs to start in its basin of attraction, and the 
rate of approach will asymptotically depend on the local linearization characteristics of our fixed point formulation. It 
is therefore desirable to develop approaches that do not rely on the stability of the fixed point, and the Newton method 
is the obvious choice. Lack of an explicit macroscopic equation means that the Jacobian involved in Newton iterations 
is not explicitly available. In principle one could estimate the Jacobian using finite differences; yet, especially for large 
scale and noisy problems, such an estimation will be both prone to error and very costly. Krylov-subspace methods, 
such as GMRES, have been devised for the iterative solution of linear equations; they are based on the evaluation 
of matrix- vector products of the system Jacobian with a sequence of algorithmically determined vectors. When the 
Jacobian is not explicitly available, matrix vector products can be estimated in a matrix -fre e fashion from nearby 
function evaluations - this is the basis for matrix- free Newton-Krylov-GMRES algorithms [jj]]- Such algorithms are 
naturally suited to the fixed point problems arising in our equation-free rcnormalization scheme - lifting from nearby 
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two-point correlation functions, evolving through the KPZ dynamics, and rescaling the resulting averaged two-point 
correlations results in an estimate of the action of the Jacobian of our fixed point problem. This type of iteration is 
often particularly well suited for problems with a separation of time scales ■ 

We used Newton-Krylov GMRES iteration to find the shape of the fixed point of the KPZ renormalization problem; 
the operating parameters were L = 400, A = 5, and a short burst of time simulation, At = 5, was used to calculate the 
evolved shape. Cubic spline interpolation was again used in the construction of the function G(k). All the data were 
averaged over 96000 replica initializations consistent with this G(k). The initial guess was obtained by starting from 
a flat interface and evolving for time of to = 5. As shown in Fig. 4, it took roughly three iterations for the plotted 
norm of the residual to decrease by one order of magnititude. The resulting fixed point is visually indistinguishable 
from the one arrived at by direct substitution, and the same scaling exponents, within the fluctuation bounds, are 
recovered. The Newton-based process is not computationally efficient in problems such as this, where the self-similar 
solution is strongly attracting; it would become advantageous, however, in cases where the self-similar solution is very 
slowly attracting, or simply unstable, when direct iteration will not converge at all. 




FIG. 5: GMRES iterations to find the fixed shape (L = 400, A = 5, A = 5, N = 96000) The inset shows the residual R as a 
function of iteration number. 



V. DISCUSSION AND CONCLUSIONS 



We presented a computational approach to finding self-similar solutions for the statistics of a KPZ-type stochastic 
evolution equation. This was accomplished without an explicit, closed form of an equation governing the dynamics of 
these statistics (in particular, the two-point correlation function); we estimated this unavailable equation on demand 
using short bursts of appropriately initialized direct simulations. The approach was successful in reproducing the 
known scaling behavior of the model. One of the important features of the computation was the detection of the 
signature of Edwards- Wilkinson scaling in the data when the working scales were not chosen large enough (du e to 
asymptotic self-similarity) . In our current implementation we used two local conditions ( "templates" , |16L Il9j ) to 
implement the two rescalings, those of the stretching of the wavenumber (or reshrinking of the length scale) and 
the scaling down of Gk (equivalently the field variable amplitude Kk)- It would be preferable to employ non-local 
conditions, making the computation more robust to local fluctuations through averaging. 

Perhaps the most important assumption of the equation-free approach is that an equation exists and closes at a 
chosen level of description; here we assumed that the appropriate level was the two-point correlation function, and 
that higher order correlations either do not affect this evolution or become quickly slaved to it. It is known |2jj that in 
one spatial dimension the exponents are insensitive to variations in the third and higher order correlators; this appears 
not to be true in higher dimensions. In our computational experiments we found that third order correlations did not 
become quickly slaved to two-point correlations (over times comparable to the evolution of the two-point correlation 
itself); since we worked in 1+1 dimension, this is consistent with the above observation. In higher dimensions, 
it is not clear that an equation does indeed close in terms of only the two-point correlation function; testing this 
hypothesis would be an important first task to pursue with our approach. Computational approaches to initializing 
"fast" variables consistently with slow ones (alternatively, on a manifold parametrized by the slow ones) have long 
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been known in computational chemistry ( 21, 22]), and their use in an equation-free context is discussed, for example, 
in H. 

In this paper we used a KPZ-type SDE as our "inner" , fine scale solver. The procedure is identical if the SDE 
solver is substituted by, for example, a kinetic deposition model; the computation of coarse self-similar solutions for 
such models is underway. For some of these, such as ballistic aggregation, we do not anticipate any difficulties; for 
other, more highly constrained models, e.g. the restricted SOS model y| the lifting operation should be nontrivial. 
In the case of stable self-similar solutions additional equation-free techniques, such as coarse projective integration 
[24II2H can be used to accelerate the computation of self-similar dynamics. 

An important test of the approach will be its ability to compute and characterize the unstable fixed point that is 
known to exist in d > 3 dimensions; this is a case where the Newton-GMRES procedure would be crucial. 
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